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Abstract 


Two  methods  for  computing  local  mass  flux  for  a  continuous  Galerkin  finite  element  formulation  of  the  Generalized 
Wave  Continuity  Equation  (GWCE)  are  derived  and  a  third  method  is  discussed  in  light  of  the  first  two.  The  GWCE 
is  shown  to  not  conserve  mass  locally,  while  it  can  be  shown  to  conserve  a  certain  quantity  locally.  The  two  derived 
methods  are  demonstrated  for  a  realistic  tidal  flow  problem  in  the  Bight  of  Abaco,  Bahamas. 

Published  by  Elsevier  B.V. 

Keywords:  Galerkin  finite  element;  GWCE;  Local  conservation;  Local  mass  flux 


1.  Introduction 

The  need  to  preserve  conservation  within  the  discrete  representation  of  a  governing  conservation  law  is 
not  in  dispute.  There  are  however,  several  ways  by  which  the  conservative  properties  of  a  particular  algo¬ 
rithm  can  be  evaluated.  For  more  than  a  decade,  the  locally  conservative  properties  of  continuous  Galerkin 
(CG)  finite  element  formulations  of  the  shallow  water  equations  have  been  in  question,  see  [16,12,2].  In  par¬ 
ticular  mass  imbalances  were  first  observed  independently  in  1992  by  Kolar  and  Westerink  [9]  in  transient 
solutions  obtained  using  the  Generalized  Wave  Continuity  Equation  (GWCE)  form  of  the  primitive  mass 
conservation  law.  These  collective  experiences  have  led  to  the  conclusion  that  CG-based  finite  element  solu¬ 
tions  are  not  locally  conservative.  The  mass  error  used  to  evaluate  local  conservation  in  these  applications 
to  date  has  typically  been  computed  element-wise  using  a  volumetric  approach  by  which  the  discrete  form 
of  the  primitive  conservation  law  is  formulated  using  physical  arguments,  see  [8]  for  example.  While 
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conceptually  satisfying,  this  approach  is  inconsistent  with  the  finite  element  basis  and  equations  used  to 
solve  the  GWCE  itself  as  brought  to  light  by  Berger  and  Howington  [3]  and  Hughes  et  al.  [6].  In  their  work, 
Berger  and  Howington  [3]  show  how  to  define  partial  fluxes  in  such  a  way  as  to  be  consistent  with  the  local 
CG  conservation  statement.  They  also  argue  that  with  this  definition  of  the  partial  flux,  the  CG  method  can 
be  seen  as  locally  conservative  when  applied  to  conservation  laws.  In  recent  work  by  Hughes  et  al.  [6],  they 
derive  an  appropriate  post-processing  technique  that  through  expansion  of  the  basis  of  the  solution,  results 
in  a  CG  finite  element  solution  that  is  locally  conservative,  i.e.  on  each  element  the  conservation  statement 
being  approximated  is  satisfied  exactly  by  the  discrete  solution.  Hughes  et  al.  [6]  demonstrate  that  the  CG 
method  is  not  only  globally  conservative  but  is  also  locally  conservative  when  using  a  strictly  consistent 
means  for  evaluating  conservation.  To  be  clear,  when  checking  mass  balance,  we  say  that  a  method  is  con¬ 
sistent,  if  it  enforces  mass  balance  in  the  same  manner  that  is  prescribed  by  the  discrete  equations  which  are 
solved,  this  includes  using  the  same  set  of  equations  and  the  same  basis  functions  as  the  solution  itself  uses. 

One  important  distinction  regarding  the  analyses  of  both  Hughes  et  al.  [6]  and  Berger  and  Howington  [3] 
is  that  they  are  relevant  to  the  primitive  form  of  the  continuity  equation.  Of  particular  interest  here  is  local 
conservation  with  respect  to  a  GWCE-based  shallow  water  equation  model.  So,  in  light  of  the  work  of 
Hughes  et  al.  [6]  and  Berger  and  Howington  [3],  we  endeavor  to  examine  the  local  conservation  properties 
of  the  GWCE  by  constructing  fluxes  that  are  consistent  with  the  CG  statement  for  the  GWCE,  which  is  not 
an  explicit  conservation  law.  In  doing  so  we  show  that  the  GWCE  is  not  locally  mass  conservative, 
although  it  is  locally  conservative  with  respect  to  certain  quantities. 

The  remainder  of  the  paper  is  organized  as  follows:  In  Section  2  we  briefly  go  through  the  derivation  of 
the  GWCE  and  cast  it  into  an  element  level  statement  of  the  problem.  In  Section  3  we  derive  and  discuss 
three  ways  of  computing  local  mass  fluxes.  In  Section  4  we  give  a  numerical  example  of  our  local  mass 
fluxes  and  the  resulting  local  mass  balance  calculations.  We  conclude  the  paper  with  a  summary  of  our  find¬ 
ings  and  acknowledgments. 


2.  Derivation  of  the  GWCE  for  ADCIRC 

We  wish  to  derive  a  consistent  method  of  computing  mass  flux  for  the  two-dimensional  coastal  circula¬ 
tion  model  known  as  ADCIRC,  see  [11,13].  ADCIRC  solves  the  generalized  wave  continuity  equation 
instead  of  the  primitive  form  of  the  vertically-integrated  continuity  equation  which  is  given  as 

^+I((/h)+I(w)=o'  "> 

where  (7,  V  —  jj  w,  vdz  is  the  depth-averaged  velocities  in  the  x,  y  directions;  w,  v  is  the  vertically-varying 
velocities  in  the  x,  y  directions;  H  -  (  +  h  is  the  total  water  column  thickness;  h  is  the  bathymetric  depth;  £  is 
the  free  surface  departure  from  the  geoid. 

ADCIRC  [1 1]  and  other  CG  finite  element  coastal  circulation  models,  such  as  Lynch  et  al.  [14],  use  the 
GWCE  instead  of  Eq.  ( 1 )  as  a  means  of  controlling  spurious  2Ax  oscillations  that  can  result  from  identical 
basis  functions  being  used  for  both  elevation  and  velocity  approximations,  see  [14,10].  One  of  the  tradeoffs 
is  that  local  mass  conservation  is  no  longer  explicitly  enforced  by  the  Galerkin  finite  element  statement. 
This  does  not  present  a  problem  in  all  situations  of  interest  to  the  modeler,  but  when  one  wants  to  couple 
the  results  from  the  circulation  model  to  a  transport  model,  local  mass  conservation  is  a  must.  Also,  for  3D 
computations,  incorrect  vertical  velocity  profiles  can  be  attributed  to  the  lack  of  local  mass  conservation  in 
the  horizontal  current  fields,  see  Luettich  et  al.  [12]. 

What  follows  is  a  summary  of  the  theoretical  formulation  of  the  GWCE,  originally  developed  by  Kinn- 
mark  [7],  as  formulated  and  solved  in  the  2D/3D  version  of  ADCIRC.  The  details  of  the  derivation  follow 
closely  those  presented  by  Luettich  and  Westerink  [13]. 
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To  start  we  assume  that  we  have  a  horizontal  computational  domain  Q  C  U2  and  its  boundary(ies)  f  > 
from  which  we  create  a  two-dimensional  unstructured  finite  element  mesh  consisting  of  N  nodes  and  M 
non-overlapping  triangular  elements  An  such  that  no  element  crosses  the  physical  boundary(ies)  T  and 
\J^_  |d„  =  (2.  Additionally  at  each  node  we  are  given  a  bathymetric  depth  h  that  varies  linearly  between 
adjacent  nodes.  In  order  to  derive  the  GWCE  we  begin  by  first  differentiating  Eq.  (1)  with  respect  to  time 
and  assuming  a  time  invariant  bathymetric  depth  (i.e.  ™  =  |£),  which  leads  to 

0  +  ^(W»  +  4(W)=O-  (2) 

Next  we  multiply  Eq.  (1)  by  the  parameter  t0  which  may  be  variable  in  space,  and  add  the  result  to  Eq. 
(2)  to  get 


We  use  the  chain  rule  and  rearrange  terms  to  get 

02c  0C  c )JV  0TO  0TO 

+  +  +  vh^  =  o, 

or  of  ox  oy  0x  oy 


(3) 


(4) 


where 


J„  =-(UH)  +  t0UH 

a  a 

=ir+'”a 

=  H^  +  U~+t„UH, 
0/  0/ 


(5a) 

(5b) 

(5c) 


9a 
a  i 


+  r«a 


3K  8C 

=  /y  a7  +  Ka/  +  r°w’ 


(6a) 

(6b) 

(6c) 


and  0A1g,.  =  UHy  VH  is  the  x, ^-directed  flux  per  unit  width. 

A  weighted  residual  method  is  applied  to  Eq.  (4)  by  multiplying  each  term  by  a  weighting  function  </>y,  to 
be  defined  later,  and  integrating  over  the  horizontal  computational  domain  Q 


(7) 

We  then  integrate  by  parts  the  terms  involving  J  x  and  jv  and  use  Eqs.  (5b)  and  (6b)  to  get  a  weak  form  of 
Eq.  (7), 

,  9r0 


j  d±\  Jj  a*A 

3/2  ’  T7„  ‘  \'"  0/  ’  "Vo  \  dx  )u  \  y'  dy  ja 


S’4„+('4'^/0 


j,  (a^+Toa)a9^'  =  o- 


(8) 


where  QN  =  [QXJ  Qy]  ■  N  is  the  outward  flux  per  unit  width  normal  to  the  boundary  r. 
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ADCIRC  in  its  present  formulation  makes  use  of  the  non-conservative  form  of  the  vertically-integrated 
momentum  equations  given  as 


dU  dU  dU  n 

- h  U  +  V  ~z - fV 

dt  9jc  dy 


dV  SK  dV  „ 

-+u  —  +v—-JV  =  -g- 
9  /  dx  dy 


r  p s 

0 

t,  + - at/ 

_ 

sa>  J 

7 

(,  H - a/7 

gPa 


0V 


HPo 


T,c 


T/u- 

A/r 

Dx 

By 

#Po 

+  77 

H 

”77 

% 

A/,. 

Dy 

By 

WPo 

+  7T 

~7 7~ 

"77' 

(9a) 

(9b) 


where /is  a  Coriolis  parameter,  g  is  the  acceleration  due  to  gravity,  p0  is  the  reference  density  of  water,  a  is 
the  effective  Earth  elasticity  factor,  r]  is  the  Newtonian  equilibrium  tidal  potential,  Ps  is  the  atmospheric 
pressure  at  the  sea  surface,  iyv,  tiT  are  imposed  surface  stresses,  i/7V,  x hy  are  bottom  stress  components, 
Mx,  Mv  are  vertically-integrated  lateral-stress  gradients,  Z)x,  Dv  are  momentum  dispersion,  5V,  Bv  are 
vertically-integrated  baroclinic-pressure  gradients. 

Now,  if  we  substitute  Eqs.  (9)  into  Eqs.  (5c)  and  (6c)  and  isolate  the  linear  free  surface  gravity  wave 
terms  for  computational  reasons,  see  for  example  Haidvogel  [4],  we  get 


9C 

A' 

Jv  —Jv  —  gh^r , 
dy 


(10a) 

(10b) 


where 


dU 


dU 


g 


J'=-e*t+-Q'i»+iQ'2t-8H- 

+  M x  —  Dx  —  Bx  4-  U  —  ■+■  to{9  , 

0/ 


dv  dy  g  a l~ 

/  =  -&*-  +  ~Qvl-  +  fQs  A  ~gH' 

ox  o  y  2  oy 


A 

.gPa 


OLY] 


dx 


A 

.gPa 


CLtf 


dy 


Pa  Pa 


T  Vi- 

Pa  Pa 


(1  la) 


+  MV-  Dv-  Bv+  y^  +  x0Qr. 

Ot 


(iib) 


Substituting  Eqs.  (10)  into  Eq.  (8)  and  rearranging  terms,  we  arrive  at  the  weighted  residual  weak  form 
of  the  GWCE  that  is  solved  by  ADCIRC, 


di~  ’ 

/ 


\a 


Q 

9t() 

"  aSc 


+  (  T0 


8C 

< 

9/ 


+  (gA 


ac 

0x’  0X 


+ 


0£  00A 


.0, 


0y’  0y 

+  /  (AA0r  =  °' 


A, 


d<t>j 

dx 


Jy 


dA\ 

'  9vA 


(12) 


We  recast  Eq.  (12)  from  a  global  statement  to  a  local,  element-wise,  statement,  by  considering  node  j 
along  with  all  nodes  adjacent  to  it  and  the  horizontal  weighting  function  </>;.  The  weighting  function  </>, 
is  defined  to  be  equal  to  one  at  node  j  and  equal  to  zero  at  all  other  nodes  such  that  varies  linearly 
between  node  j  and  the  nodes  adjacent  to  it,  see  Fig.  1 .  For  future  reference  we  define  a  finite  element  basis 
space  SPN  to  be, 

=  span {0,.y  =  1  :  N}  c  JT'(£2), 


(13) 
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Fig.  1.  Node  j  with  adjacent  nodes  and  the  associated  weight  function  (f)j. 


where  J^](Q)  is  the  usual  Sobolev  space  consisting  of  all  functions  such  that  the  function  and  its  first  deriv¬ 
atives  are  both  square  integrable  over  the  domain  Q. 

Two  other  important  properties  of  the  weighting  function  are  that  for  any  element  with  local  node  num¬ 
bers  i  =  l,  2,  3  we  have  that  the  weight  functions  over  that  element  sum  to  unity  and  the  sum  of  the  spatial 
derivatives  of  the  weight  function  over  the  element  is  zero, 

<t>i  = 1  > 

r  I 

y-  90/  _  Q  _  90, 

dx  0 y 

i  i  t= i  ^ 


(14a) 

(14b) 


At  this  time  we  define  some  additional  notations  associated  with  node  j,  A„  is  the  area  of  element  n,  NEj 
is  the  number  of  elements  containing  node  y,  AJn  is  the  local  element  number  n  which  contains  node 
j,  Qj  =  is  the  union  of  all  the  elements  containing  node  y,  ry  is  the  boundary  of  Qn  Pn  is  the 

total  boundary  of  local  element  number  n ,  ifl  is  a  two  edge  inner-element  partial  boundary  of  element 

Aj. 

n 

Using  this  new  notation  and  the  definition  of  0/  we  see  that  the  integrations  performed  in  Eq.  (12)  over 
the  entire  domain  Q  can  be  reduced  to  integrations  over  the  patch  Qj,  which  corresponds  precisely  with  the 
region  over  which  0,  ^  0, 


a  , 
Qx  dx  ’  ^ 


Q, 


Q, 


+ 

djo 

dy 


'd£  d(f)j 
1 5,v '  dx  / 

AT 


dA\ 

dy] 


Q, 


(!>Q„ 

l“aT 


T  To Qn  )  0y9r j  —  0. 


(15) 


We  note  that 

ft  Ci^  +  To0/v)0/9Ty  =  0  whenever  node  y  does  not  lie  on  the  physical  boundary  f,  thus  Eq.  (15)  is 
equivalent  to  Eq.  (12).  Also  due  to  cancellation 
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NE. 


( +  t»2,v  )  4 fir,  =  Y 


6/ 


+  *<>(?, v  <M/' 


(16) 


So  we  rewrite  Eq.  (15)  using  elemental  notation  as 


f  M 


NEj 

E 


i?A')K + {A*),, + AiA)K + (AW, 
AW)  A  A)  AAA,,. 
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which  is  the  element  level  spatial  discretization  of  the  GWCE  with  respect  to  the  weight  function  4>j. 


(17) 


3.  In  search  of  a  consistent  and  conservative  flux 


In  this  section  we  formulate  two  ways  of  computing  local  mass  fluxes  and  the  resulting  computation  of 
mass  error.  We  then  discuss  a  third  approach  in  light  of  the  other  two. 


3.L  Finite  volume  flux 


We  begin  with  a  simple  and  straightforward  way  of  computing  local  mass  fluxes  and  the  resulting  com¬ 
putation  of  mass  error.  This  method  is  commonly  used,  see  for  example  [8,5,2],  easy  to  compute  and  gives  a 
physically  realistic  mass  flux,  in  that  it  is  defined  normal  to  an  element’s  edge.  The  local  mass  flux  defined  in 
this  section  is  akin  to  one  used  by  Kolar  et  al.  [8]  except  we  do  not  integrate  in  time. 

We  introduce  a  new  space  SF  C  jS?2(I2),  such  that  Sf  —  span{w„  :  n  —  1  :  M }  where  w„(x,y)  is  equal  to 
one  wherever  {x,y)  e  A„  and  zero  everywhere  else.  We  begin  with  a  weak  statement  of  the  primitive  form  of 
the  continuity  equation  (1)  over  an  arbitrary  element 


SC 
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where  wn(x,y)\A  =  1.  Next  we  integrate  the  second  and  third  terms  by  parts  and  rearrange  to  get, 


L 


QA,w„dr „  = 


An 


Thus,  we  define  the  first  of  our  local  mass  fluxes,  what  we  will  call  a  finite  volume  flux,  as 
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QNwndr„. 


(18) 


(19) 


(20) 


We  can  evaluate  this  mass  flux  directly  over  any  edge  of  any  element  using  the  computed  finite  element 
solution.  This  mass  flux  is  uniquely  defined  by  the  nodal  values  at  the  endpoints  of  an  element  edge,  which 
means  that  the  edge  flux  from  neighboring  elements  will  differ  only  in  sign,  due  to  the  outward  normals. 
However,  this  method  is  not  consistent  with  the  finite  element  statement  solved  in  ADCIRC  in  that  it  is 
derived  from  Eq.  (1)  rather  than  Eq.  (17)  which  is  solved  by  ADCIRC,  furthermore,  the  weight  function, 
w„  e  is  not  part  of  the  original  basis  space,  .  As  a  result  local  mass  balance  in  the  context  of  this  flux 
is  not  guaranteed,  by  this  we  mean  that  if  we  evaluate  Eq.  (19)  directly  from  the  GWCE-based  computed 
solution,  we  will  not  in  general  obtain  equality. 
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This  method  remains  attractive  because  the  resulting  fluxes  are  physically  realistic,  it  is  easy  to  compute, 
and  as  Kolar  et  al.  [8]  note,  it  is  a  good  indicator  of  a  methods  ability  to  conserve  mass. 

3.2.  A  consistent  flux 


Since  the  finite  volume  flux  is  inconsistent  with  our  finite  element  statement,  we  also  consider  a  local 
mass  flux  that  is  defined  in  a  manner  consistent  with  the  finite  element  formulation  of  the  GWCE.  This 
approach  is  similar  to  one  proposed  by  Berger  and  Howington  [3]  for  conservation  laws  in  primitive  form. 

We  seek  to  derive  a  consistent  elemental  flux  from  the  discrete,  local  form  of  the  GWCE  given  in  Eq. 
(17).  In  a  fashion  similar  to  that  detailed  by  Berger  et  al.  [3],  we  consider  a  node  j  and  its  adjacent  neighbors 
which  form  the  patch  Qj.  We  present  in  Fig.  2  an  exploded  view  of  the  patch  Qj,  where  we  have  extracted 
element  A\  in  order  to  show  the  edge  fluxes  along  the  inner  element  boundaries,  r\.  Note  that  for  conve¬ 
nience  we  are  using  local  element  numbering. 

We  separate  out  the  terms  involving  element  A\  in  Eq.  (17)  and  get 
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Notice  that  since  =  0  on  Fj,  evaluating  the  boundary  integral  in  Eq.  (21)  over  f{  is  equivalent  to  eval- 
uating  it  over  just  f\. 

By  examining  the  first  bracketed  expression  in  Eq.  (21),  we  can  ascertain  the  inherent  definition  of  a  flux 
like  quantity  from  the  finite  element  statement,  similar  to  what  Berger  and  Howington  [3]  observed  in  their 


Fig.  2.  The  patch  O'  with,  A\  separated  to  show  edge  fluxes, 
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work.  Thus,  we  define  a  partial  numerical  “flux”  P]  and  its  time  rate  of  change  of  p\  through  the  interface 
as 


K,  +  A,= 
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ac  0</>, 
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(  Q"  9 y  ’ ,, 


(22) 


Here  we  assume  that  node ydoes  not  lie  on  the  boundary  r,  otherwise  we  must  adjust  Eq.  (22)  so  that  we 
use  the  appropriate  boundary  conditions  set  forth  in  the  model  specifications.  We  recognize  the  right  hand 
quantity  is  a  function  of  time  and  thus  treat  Eq.  (22)  as  an  ODE  that  we  numerically  solve  using  Huen’s 
method,  which  is  second  order  accurate.  For  discussion  purposes  we  present  the  analytical  solution 
such  that 
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and  c  =  Pj  (0).  For  most  applications,  the  water  surface  is  initially  at  rest,  so  that  c  -  0.  The  quantity  P\{l) 
represents  only  a  partial  mass  flux  for  element  A\ ,  namely  the  contributions  to  the  total  flux  associated  with 
node  j.  Thus  to  find  the  total  mass  flux  for  element  A\,  we  sum  the  partial  fluxes  associated  with  the  vertices 
of  that  element,  Py(t)  =  />’.  (/).  It  is  not  clear  if  this  partial  flux  is  physically  realistic,  certainly  one  can 

not  get  the  flux  across  a  single  edge  of  an  element  with  this  method. 

Notice  that  the  solution  of  (23)  is  over  the  time  interval  [0,  t]  while  we  are  interested  in  the  solution  in  the 
time  interval  To  obtain  a  time  integration  over  the  desired  time  interval  we  subtract  subsequent 

solutions,  namely, 


p'B(tk)  =  p'(h)-p\h->), 


(25) 


and  define  /^(/*)  as  a  local  mass  flux  for  element  A\  at  time  tk. 

To  complete  the  derivation  of  F'B  and  to  ensure  consistency,  we  first  must  use  the  integration  rules 
employed  by  ADC1RC  for  the  components  of  <P//),  given  in  Table  1,  whereby  *Pj{t)  becomes 
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(26) 


Finally  we  note  that  in  solving  the  GWCE,  ADCIRC  uses  a  weighted  three  step  semi-implicit  time  dis¬ 
cretization  procedure  at  s  —  I,  s,  and  s  +  1.  All  the  time  steps  are  uniform  in  size,  i.e.  At  =  ts  —  The 
time  dependent  quantities  in  Eq.  (26)  are  approximated  in  the  following  way 
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Table  1 

Relevant  integration  rules  used  in  ADCIRC 


=T*' 

(^)A 


2  °' 

n  0To  j.  \  r,  v-3  h' 


Qi 


Qy-fr'*})  =  QviE;_ito/ 


9  2c,  c;+l-2  c;  +  c;' 
s  r-  m2 

(27a) 

k  _  c  -  cr' 

3/  2A  /  ’ 

(27b) 

£•  =  «|C’  +  °bCJ  +  «jC'» 

(27c) 

Jx  |  ,  J  y  |  =Jx\l*ty\' 

(27d) 

Q^-Qy\  =  &vQ\v 

(27e) 

where,  cq,  a2,  23,  are  the  time  weighting  factors  which  are  set  by  the  user. 

In  terms  of  realizing  perfect  local  mass  conservation  with  this  consistent  flux,  it  is  worth  looking  at  the 
analytical  solution  to  the  total  flux  for  element  A\,  prior  to  time  discretization,  i.e.  Eq.  (23).  We  assume  for 
the  time  being  that  r0  is  constant  on  the  element;  this  assumption  is  not  prohibitive  and  is  presently  the 
default  case  for  historical  simulations  using  ADCIRC.  For  each  element,  we  sum  the  partial  mass  flux  asso¬ 
ciated  with  each  node  to  get 


/5l(/)  =  E<w  =  Ee' 
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1,v  (s)ds  =  e 
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Next  we  use  the  definition  of  ^j(/)  and  Eq.  (14a)  and  (14b)  to  get 


Now  we  apply  Leibniz’s  rule  to  switch  the  order  of  integration  and  differentiation  along  with  the  assump¬ 
tion  that  on  the  current  element,  r0  is  constant,  to  obtain 


We  recognize  that  the  integrand  is  a  result  of  a  product  rule  differentiation,  so  upon  rearrangement  we 
get 
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(31) 
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By  the  Fundamental  Theorem  of  Calculus,  Eq.  (31)  collapses  to  a  statement  of  local  mass  conservation, 
namely 

p,{,)= -(%■'),{  (32) 

Of  course,  when  we  do  discretize  Eq.  (23)  in  time  whereby  we  can  no  longer  exactly  perform  the  operations 
done  to  get  Eqs.  (30)-(32),  we  should  not  expect  perfect  local  mass  conservation.  This  is  mainly  due  to  the 
presence  of  numerical  errors  in  the  approximation  of  the  second  derivative  of  elevation  with  respect  to  time. 

In  Section  3.1  we  developed  a  flux  that  is  physically  realistic,  yet  is  not  consistent  with  the  finite  element 
formulation  of  our  problem  and  hence  does  not  guarantee  local  mass  conservation.  In  Section  3.2  we  devel¬ 
oped  a  flux  that  is  consistent  with  the  finite  element  formulation,  but  may  not  be  physically  realistic  and  still 
does  not  result  in  local  mass  conservation  due  to  the  required  time  discretizations.  This  leads  us  to  examine 
another  option  for  determining  a  consistent  local  mass  flux  for  the  GWCE. 


3.3.  A  post-processed  flux 

The  final  flux,  a  post-processed  flux,  will  only  be  discussed  in  the  context  of  the  GWCE  and  follows  the 
ideas  set  forth  in  the  work  of  Hughes  et  al.  [6].  In  their  work,  Hughes  et  al.  [6]  show  another  method  for 
defining  local  fluxes  in  such  a  way  as  to  be  consistent  with  the  finite  element  state,  see  also  Oshima  et  al. 
[17].  As  a  result  they  are  able  to  show  that  the  continuous  Galerkin  finite  element  method  when  applied 
to  conservation  laws  like  Eq.  (1),  is  locally  conservative.  This  is  accomplished  through  a  post-processing 
of  the  original  finite  element  statement.  Their  approach  dispels  a  common  and  long  held  notion  that  the 
CG  method  is  only  globally  conservative. 

Their  basic  idea  is  to  solve  the  original  finite  element  statement  first,  then  define  an  auxiliary  flux  for  the 
global  boundary,  whereby  they  solve  a  modified  finite  element  problem  for  this  auxiliary  flux.  Once  this  is 
complete  they  solve  modified  local  finite  element  problems  for  element  level  fluxes.  The  modifications  con¬ 
sists  of  adding  new  boundary  terms  for  the  areas  of  interest,  elements  in  our  case.  Furthermore,  the  original 
basis  space  £fN  is  enlarged  by  including  all  the  weight  functions  wn  in  Sf  ,  when  solving  for  element  level 
fluxes.  Using  this  procedure,  Hughes  et  al.  [6]  show  that  the  resulting  fluxes  are  locally  and  globally  con¬ 
servative.  One  consequence  of  solving  the  additional  problem  for  the  fluxes,  they  point  out,  citing  Babuska 
and  Miller  [1],  is  that  these  new  fluxes  have  many  desirable  qualities  about  them,  including  superconver¬ 
gence  properties.  Their  work  falls  into  a  broader  category  of  methods  called  “goal-oriented  error  estima¬ 
tion”  or  “quantities  of  interest”  see  for  example  [19,18,15],  where  estimates  are  obtained  for  integral 
functionals,  such  as  mass  flux  in  our  case  and  lift  or  drag  in  aeronautical  applications  to  list  just  a  few. 

By  applying  the  ideas  detailed  in  Hughes  et  al.  [6]  to  the  GWCE  in  an  effort  to  preserve  conservation  and 
consistency,  we  need  to  solve  a  system  of  equations  for  not  only  the  mass  flux  but  also  the  time  rate  of 
change  of  the  mass  flux,  e.g.  quantities  similar  to  those  on  the  left  hand  side  of  Eq.  (22).  Both  of  these  quan¬ 
tities  taken  together  are  necessary  to  achieve  local  conservation  and  maintain  consistency  with  the  GWCE- 
based  finite  element  formulation.  But  what  we  desire  is  only  the  local  mass  flux.  So  a  similar  situation 
presents  itself  to  the  one  we  had  in  Section  3.2  where  an  ODE  was  solved  in  order  to  extract  the  mass  flux 
quantity.  The  need  to  discretize  in  time  would  again  result  in  a  lack  of  local  mass  conservation. 

As  a  word  of  caution,  Hughes  et  al.  [6]  point  out,  that  their  post-processing  technique  is  basically  an  L2 
projection  which  can  result  in  spurious  oscillations  near  discontinuities.  If  spurious  oscillations  are  reintro¬ 
duced  using  such  a  post-processing  technique,  then  we  have  defeated  the  purpose  of  using  the  GWCE  in  the 
first  place. 

It  is  possible  to  use  the  ideas  of  Hughes  et  al.  [6]  to  setup  the  post-processing  procedure  to  solve  for  just 
the  mass  flux,  i.e.  pretend  we  are  solving  Eq.  (1).  This  approach  results  in  a  locally  conservative  mass  flux, 
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but  is  not  consistent  with  the  original  equation  solution  from  the  GWCE  and  very  well  may  reintroduce 
spurious  oscillations.  Therefore,  even  considering  the  method  of  Hughes  et  al.  [6],  we  are  not  able  to  satisfy 
both  consistency  and  local  mass  conservation  for  a  GWCE-based  model. 

While  the  methodology  laid  down  by  Hughes  et  al.  [6]  for  finding  locally  conservative  fluxes  is  valuable, 
the  resulting  flux  is  relatively  expensive  to  compute  and  is  an  element  averaged  quantity  which  is  not  useful 
for  transport  purposes.  Furthermore,  the  lack  of  local  mass  conservation  caused  by  the  need  to  solve  an 
ODE  for  the  flux  and  the  possibility  of  reintroducing  spurious  oscillations  through  Eq.  (I)  are  counterpro¬ 
ductive.  Thus,  for  these  reasons  we  do  not  derive  the  fiux  as  proposed  by  Hughes  et  al.  [6]  for  the  GWCE, 
nor  give  any  numerical  examples  of  it. 


l  2  3  4  5  6  7  1234567” 

x  (m)  x  jq4  x  (m)  x  jq* 


(a)  (b) 

Fig.  3.  The  (a)  mesh  and  (b)  bathymetry  in  meters,  for  the  Bight  of  Abaco,  Bahamas  with  selected  stations  marked  with  a  <g>. 


Fig.  4.  For  the  Bight  of  Abaco,  Bahamas:  (a)  Computed  vorticity  of  the  flow  (I/s)  at  day  12.0  of  the  simulation  with  time-series 
stations  marked  with  a  ®  and  (b)  the  lime-series  of  vorticity  at  each  station  from  day  1 0.0  to  12.0. 
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4.  Numerical  example 

The  finite  volume  flux,  /Vv,  and  the  consistent  flux,  PB ,  defined  in  the  previous  sections  for  the  GWCE- 
based  finite  element  method  employed  by  the  coastal  ocean  circulation  model  ADCIRC,  have  been  tested 
for  a  realistic  coastal  flow  problem.  The  results  of  this  test  are  presented  below  in  an  effort  to  support  the 
analytical  findings. 

We  present  a  realistic  application  for  tidal  circulation  in  the  Bight  of  Abaco,  Bahamas,  which  lies  be¬ 
tween  the  islands  of  Great  Abaco  to  the  east,  Little  Abaco  to  the  north,  and  Grand  Bahama  to  the  west. 
This  region  is  a  site  of  significant  nonlinear  tidal  generation  due  to  its  isolation  from  deep  offshore  waters 
(1000-2000  m)  and  its  own  very  shallow  depths.  Tides  propagate  from  open  waters  and  are  largely  modified 
by  frictional  dissipation  and  nonlinear  tidal  interactions  within  the  Bight. 


Fig.  5.  For  the  Bight  of  Abaco,  Bahamas:  Computed  (a)  elevation  (m),  (b)  magnitude  of  velocity  (m/s),  (c)  magnitude  of  the  .v 
component  of  velocity  (m/s),  and  (d)  magnitude  of  the  y  component  of  velocity  (m/s)  at  day  12.0  of  the  simulation  with  selected  stations 
marked  with  a  <g>. 
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The  basin  is  discretized  by  an  unstructured  triangular  mesh  using  1696  elements  and  926  nodes,  see  Fig. 
3a,  such  that  the  maximum  element  radius  is  1750  m  and  the  minimum  element  radius  is  454  m  (element 
radius  being  defined  as  the  radius  of  the  circle  that  circumscribes  the  element).  The  bathymetry  varies  from 
only  1 .0  to  9.0  m  as  shown  in  Fig.  3b.  Tidal  forcings  consisting  of  five  primary  constituents  ( K\ ,  O i,  M2,  S2, 
and  yV2)  are  applied  at  the  open  ocean  boundary  and  are  used  to  drive  the  )  2-day  simulation  that  includes  a 
5  day  ramp  period.  A  fixed  time  step  value  of  30  s.  is  used  along  with  a  GWCE  parameter  set  to  t0  =  0.01. 
Lateral  mixing  effects  are  not  considered. 

For  the  Bight  of  Abaco  simulation,  the  minimum  and  maximum  elevation  values  over  approximately 
four  M2  tidal  cycles  (day  10  to  day  12)  were  -0.3622  and  0.4594  m,  respectively.  The  x  component  of  velo¬ 
city  ranges  (from  —0.42 12  to  0.3594  m/s)  and  they  component  of  velocity  ranges  (from  —0.3734  to  0.3140  m/ 
s)  with  a  resulting  maximum  magnitude  of  velocity  of  0.4451  m/s.  Fig.  5  presents  a  snapshot  from  day  12.0 


Fig.  6.  For  the  Bight  of  Abaco,  Bahamas:  Computed  fluxes  (a)  P\v,  the  finite  volume  flux,  and  (b)  Pti,  the  consistent  flux,  and  their 
corresponding  mass  accumulation  terms  (c)  and  (d)  — 1)  at  day  12.0  of  the  simulation:  selected  stations  are  marked  with  a  <g>.  Note 
that  (c)  and  (d)  are  the  same  quantities  viewed  at  different  color  scales. 
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Table  2 


Invariant  mesh  and  bathymetric  information  relevant  to  each  of  the  station  locations  for  the  Bight  of  Abaco.  Bahamas 


Station  # 

Element  area  (km)2 

Node  #  / 

Depth  hj  (m) 

II*/  -  Mk  (m) 

Nodal  radius  (m) 

Edge  length  (m) 

1 

2.7223 

412 

6.00 

1.500 

2687.4 

2608.8 

440 

4.75 

1.500 

2608.8 

2522.1 

436 

6.00 

1.250 

2687.4 

2403.7 

2 

2.7083 

792 

3.75 

0.500 

2713.1 

2534.5 

813 

3.50 

0.250 

2713.1 

2713.1 

809 

3.25 

0.500 

2760. 1 

2304.9 

3 

2.7778 

317 

7.50 

0.750 

2833.3 

2669.0 

336 

7.25 

0.500 

2713.1 

2687.4 

312 

7.75 

0.500 

2833.3 

2713.1 

Fig.  7.  For  the  Bight  of  Abaco,  Bahamas:  Comparisons  of  the  finite  volume  flux  Pfv.  the  consistent  flux  PB,  and  the  mass  accumulation 
—  1)  in  mVs  for  (a)  station  I.  (b)  station  2,  and  (c)  station  3. 

of  the  computed  elevation,  magnitude  of  the  velocity,  and  values  of  the  a-  and  _y-velocity  components  over 
the  entire  mesh,  while  Table  3  list  the  specific  values  of  these  quantities  at  the  nodes  of  the  three  stations  to  be 
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discussed  in  more  detail  later.  We  also  show  the  corresponding  computed  two-dimensional  vorticity  of  the 
flow  (1/s)  in  Fig.  4  for  day  12.0. 

For  the  same  time  period  we  present  values  of  the  local  mass  fluxes  (m 3/s),  computed  as  the  finite  volume 
flux,  P\\,  from  Eq.  (20)  and  the  consistent  flux  PB  from  Eq.  (25)  along  with  their  corresponding  mass  accu¬ 
mulation  terms  l)  in  Fig.  6.  The  color  scales  for  Fig.  6a-c  are  fixed  between  -573  and  673  m3/s  in 
order  to  facilitate  a  direct  comparison  between  the  two  flux  measures.  What  is  immediately  obvious  in 
the  computed  flux  terms  is  that  no  apparent  pattern  emerges  for  the  consistent  flux  values  PB,  while  the 
finite  volume  flux  P fv  seems  to  track  the  tidal  propagation  in  and  out  of  the  basin  as  shown  in  Fig.  5.  FIow- 
ever,  neither  flux  balances  their  respective  mass  accumulation  term  —  the  consequence  being  that 

local  mass  is  not  conserved.  Notice  that  the  largest  mass  accumulation  values  in  Fig.  6d  occur  roughly 
in  the  same  regions  where  P{v  is  largest  in  Fig.  6a.  At  these  same  locations  the  vorticity  is  also  at  its  largest 
magnitude  as  seen  in  Fig.  4a.  The  correspondence  between  vorticity  and  local  mass  imbalance  with  respect 
to  P{\,  is  discussed  in  more  detail  by  Blain  and  Massey  [2]. 


Fig.  8.  For  the  Bight  of  Abaco,  Bahamas:  Comparisons  of  the  element-averaged  (a)  elevation  (m),  (b)  magnitude  of  velocity  (m/s), 
(c)  values  of  the  .v  component  of  velocity  (m/s)  and  (d)  values  of  the  y  component  of  velocity  (m/s)  for  each  of  the  three  stations. 
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Table  3 

Information  for  the  three  stations  at  day  12.0  of  the  simulation  for  the  Bight  of  Abaco,  Bahamas 


Station  # 

Node  #  / 

Ci  (m) 

a,  (m/s) 

vf  (m/s) 

1 

412 

4.9848e-02 

-2.7433e-02 

-1.2202e-01 

440 

5.0l24e-02 

— 7.2452e-02 

—  1.2054e-01 

436 

5.2335e-02 

—3. 19 19e-02 

—  1.2376e-01 

2 

792 

l.0706e-01 

-2.411 3e-02 

— 7.3052e-02 

813 

1.0905e-01 

-3.261 5e-02 

-7.2947e-02 

809 

1 . 1 022e-0 1 

— 2.4024e-02 

-6.6203e-02 

3 

317 

4.543  le-02 

3. 1354e-03 

-1.0572e-01 

336 

4.7007e-02 

—3. 1 790e-03 

- 1 . 1093e-01 

312 

4.6062e-02 

2.3 1 8 1  e-03 

—  1.0G37e-01 

To  examine  the  time  behavior  of  these  fluxes  three  stations  are  identified,  see  Fig.  3  for  their  locations, 
based  on  the  local  mass  imbalances  computed  for  each  type  of  local  mass  flux  considered.  The  relevant 
mesh  and  bathymetric  details  for  each  station  are  recorded  in  Table  2.  Station  1  is  an  element  having  a  high 
local  mass  imbalance  when  using  the  finite  volume  flux,  P^.  By  contrast,  station  2  has  a  low  mass  imbalance 
for  either  the  finite  volume,  j Pfv,  or  the  consistent  partial  flux  PB .  Lastly,  station  3  is  an  element  that  has  a 
high  mass  imbalance  when  using  the  consistent  partial  flux  PB.  The  time  series  behavior  for  each  of  the  two 
fluxes  at  each  station  is  shown  in  Fig.  7.  Corresponding  time  series  for  the  element  averaged  quantities  of 
elevation  (m),  values  of  velocity  (m/s)  and  the  values  of  the  x  and  y  components  of  velocity  (m/s)  are  shown 
in  Fig.  8.  Again,  from  Fig.  7,  we  see  that  neither  of  the  computed  fluxes  balances  the  mass  accumulation 
term  — ,  wn)  at  any  of  the  station  locations.  Both  fluxes  appear  to  have  the  same  phase  yet  slightly  dif¬ 
ferent  shapes  at  the  extrema.  Both  fluxes  are  also  out  of  phase  with  the  mass  accumulation  term  —  (|jj ,  w„), 
especially  at  stations  l  and  3  where  the  flux  values  clearly  lead  the  mass  accumulation.  The  larger  values  of 
the  finite  volume  flux,  Piv  at  station  1,  are  physically  realistic  based  on  a  comparison  of  the  velocity  mag¬ 
nitudes  at  each  of  the  three  stations.  Station  1  has  a  velocity  magnitude  that  is  nearly  one  and  half  times  the 
magnitude  of  the  velocity  at  the  other  two  stations.  Station  1  is  also  located  within  a  region  that  experiences 
a  relatively  steep  gradient  in  bathymetry  as  depths  drop  dramatically  beyond  the  shallow  sill  at  the  entrance 
to  the  Bight  of  Abaco  (see  the  nodal  depths  for  station  1  in  Table  2).  Such  steep  bathymetric  gradients  are 
known  to  result  in  larger  mass  imbalances,  see  Oliveira  et  al.  [16]. 

This  numerical  example  illustrates  that  the  finite  volume  flux,  Pfv,  is  physically  realistic  even  though  it  is 
not  consistent  with  the  finite  element  formulation.  The  finite  volume  flux  is  also  a  good  indicator  of  where 
local  mass  imbalances  are  likely  to  occur.  Additionally,  we  have  shown  that  the  consistent  flux,  PB,  is  not 
physically  based  and  does  not  always  predict  regions  where  local  mass  imbalances  would  occur.  We  feel 
that  for  the  GWCE,  the  finite  volume  flux,  FYV,  should  be  retained  as  an  analysis  tool. 


5.  Conclusions 

We  examined  three  ways  of  computing  local  mass  flux  and  the  resulting  mass  balance  for  the  coastal 
ocean  circulation  model,  ADCIRC,  which  uses  the  GWCE  in  order  to  control  2Ax  spurious  waves.  We 
have  shown  that  none  of  the  fluxes  discussed  satisfy  all  three  desired  traits  of  being,  physically  realistic, 
locally  mass  conservative,  and  consistent  with  the  finite  element  statement  being  solved.  We  have  shown 
that  the  most  straight  forward  way  of  computing  fluxes,  the  finite  volume  flux,  Ffv,  gives  physically  realistic 
results.  We  hope  that  this  paper  has  added  clarity  to  the  concept  of  local  conservation  with  respect  to  con¬ 
tinuous  Galerkin  finite  element  statements. 
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